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ABSTRACT 

We report on an X-ray observation of the 166 Myr old radio pulsar J0108-1431 with XMM-Newton. 
The X-ray spectrum can be described by a power-law model with a relatively steep photon index 
F w 3 or by a combination of thermal and non-thermal components, e.g., a power-law component 
with fixed photon index F — 2 plus a blackbody component with a temperature of kT = 0.11 keV. 
The two-component model appears more reasonable considering different estimates for the hydrogen 
column density Njj- The non-thermal X-ray efficiency in the single power- law model is ^f^iokcV = 
Li-iok e y / E ~ 0.003, higher than in most other X-ray detected pulsars. In the case of the combined 
model, the non-thermal and thermal X-ray efficiencies are even higher, f/fj^okcV ~ wc ~ 0.006. We 
detected X-ray pulsations at the radio period of P ~ 0.808 s with significance of w 7a. The pulse shape 
in the folded X-ray lightcurve (0.15-2 keV) is asymmetric, with statistically significant contributions 
from up to 5 leading harmonics. Pulse profiles at two different energy ranges differ slightly: the profile 
is asymmetric at low energies, 0.15-1 keV, while at higher energies, 1-2 keV, it has a nearly sinusodial 
shape. The radio pulse peak leads the 0.15-2 keV X-ray pulse peak by A0 = 0.06 ± 0.03. 
Subject headings: pulsars: individual (PSR J0108-1431) — stars: neutron — X-rays: stars 



1. INTRODUCTION 

Radiation of rotation powered pulsars (RPPs) is pow- 
ered by their rotational energy loss. Over 1700 of RPPs 
have been detected in the radio and over 100 in X-rays. 
Their period and period derivative measurements are 
used to calculate the spin-down power, E, the surface 
magnetic field B, and the characteristic spin-down age, 
r. The X-ray spectra of RPPs exhibit some common 
features that evolve with spin-down age. 

For very young pulsars (r < lOkyr), the magne- 
tospheric emission is observed to be dominant, bury- 
ing the component of bulk surface thermal emission in 
most cases. Middle-aged pulsars, with spin-down ages 
10 4 < r < 10 6 yr, tend to have a significant contribution 
from the surface thermal emission (kT ~ 0.03 — 0.3 keV) 
in the soft X-rays, eventually dominating the X-ray spec- 
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trum up to ~ 1 — 2 keV. The high-energy part of the X- 
ray spectrum contains a contribution from non-thermal, 
magnetospheric emission and may contain a component 
of thermal, polar-cap emission. The non-thermal compo- 
nent in the spectra of these young pulsars is usually fit 
with a power-law (PL) with photon index 1 < T < 2 (see, 
e.g., Figure 4 in Li et al. 2008). The fraction of the spin- 
down power emitted in the X-rays (the X-ray efficiency 

?7x = Lx/ E) has values in the range of 10~ 5 < r;x < 10~ 3 
(see e.g., Figure 5 in Kargaltsev & Pavlov 2008). The 
X-ray pulse profiles of young pulsars show relatively 
sharp pulses for the non-thermal emission and smooth, 
broader pulses for the thermal emission. Asymmetric 
pulse shapes have been reported for middle-aged pulsars. 
For instance, Geminga, PSRB0656+14 and PSRB1055- 
52 (r ~ a few 100 kyr) and PSR J0538+2817 (r = 30 
kyr) show such asymmetric pulse profiles (e.g., De Luca 
et al. 2005; Zavlin & Pavlov 2004b; McGowan et al. 2003; 
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Pavlov et al. 2002). 

Old ( r Si 10 6 Y r ) R-PPs have diminished spin-down 
powers (E < 10 34 ergs _1 ) and X-ray luminosities. This 
restricts the distance up to which such old pulsars can 
be detected and necessitates longer observations. There 
have been only a dozen of these old RPPs detected in 
the X-rays (see, e.g., Posselt et al. 2012 and references 
therein). Absorbed PL model fits for these pulsars 
yield r ~ 2 — 4. Hence, their spectra are softer than 
the non-thermal spectra of the younger pulsar popula- 
tion. The inferred absorbing Hydrogen column density 
values, Nh, are usually larger than expected from the 
respective total Galactic Hydrogen column densities 
and/or estimates from pulsar dispersion measures. The 
non-thermal X-ray efficiencies in the 1-10 keV band 
show significant scatter, 771-10 keV ~ 10~ 4 — 10~ 2 , but 
are on average higher than the corresponding values 
calculated for younger pulsars (Kargaltsev et al. 2006; 
Zharikov et al. 2006). 

Old neutron stars are too cold to have significant 
thermal X-ray emission from the bulk stellar surface 
(Yakovlev & Pethick 2004). However, a possible thermal 
X-ray contribution may come from polar caps heated 
by infalling accelerated particles (Harding & Muslimov 
2001, 2002). An additional thermal component can 
explain the X-ray spectra of some of these old pulsars. 
If fitted with a simple blackbody (BB) model, these 
components have BB temperatures of 0.1-0.3 keV and 
projected areas smaller than the conventional polar-cap 
area (A pc = 2tt 2 R^ s /(cP) ~ 10 5 m 2 for NS radius 
i?NS ~ 10 km and period, P ~ Is). Thermal emission 
from the surface layers of a neutron star atmosphere can 
differ substantially from BB emission (e.g., Pavlov et al. 
1995). Assuming H or He atmosphere emission instead 
of the simple BB results in larger polar-cap area sizes 
and temperatures about a factor 2 lower than the BB 
temperature (see, e.g., Zavlin & Pavlov 2004a for the 
atmosphere model analysis of the old PSR B0950+08). 

Pulsar J0108-1431 was discovered by Tauris et al. 
(1994). It has a period P = 0.808 s and a period 
derivative P = 7.7 x 10 -17 s s _1 , which corresponds 
to r = P/2P w 166 Myr, E w 5.8 x 10 30 erg s _1 , 
and B = 2.5 x 10 11 G, as listed in the Australian 
Telescope National Facility (ATNF) Pulsar Catalogue 1 
(Manchester et al. 2005). PSR J0108— 1431 has been 
well monitored for two decades at radio wavelengths, 
and details on its radio properties have been presented 
by, e.g, Espinoza et al. (2011) and Hobbs et al. (2010, 
2004). Weltevrede & Johnston (2008) discussed the 
high linear polarization of ~ 70% at 20 cm, which is 
very unusual for a pulsar with such a low E (see their 
Figure 8). Because of suspected similarities with high-_E 
pulsars, PSR J0108— 1431 is now regularly monitored as 
part of the collaborative program with the Fermi group 
(Weltevrede et al. 2010). An optical counterpart of 
PSR J0108-1431 was suggested by Mignani et al. (2008). 

PSRJ0108-1431 has the largest r and lowest E 
among the X-ray detected pulsars. Characterization of 

1 http:/ /www. atnf.csiro.au/research/pulsar/psrcat 



its spectrum is essential to understand the final stages 
of the spectral evolution of old pulsars. It is one of 
the nearest neutron stars to Earth, with an estimated 
VLBI parallax distance of d = 240lg 24 pc (Deller et al. 
2009), which was recently revised to d = 210j^Qpc 
when accounting for the Lutz-Kelkar bias (Verbiest 
et al. 2012). Deller et al. (2009) also measured the 
pulsar proper motion of 170.0 ± 2.8 mas yr" 1 in the 
south-southeast direction. 

PSRJ0108-1431 has been detected with the Chandra 
X-ray Observatory (Pavlov et al. 2009). The spectral 
analysis, performed using the 51 counts detected in the 
0.8-5 keV energy range, showed emission consistent with 
either a single-component, soft (r ~ 2.2 — 3.4) PL ex- 
pected from the pulsar's magnetosphere or a 3 MK ther- 
mal BB emission from an apparent emitting area of ~ 50 
m 2 . Pavlov et al. (2009) suggested the presence of both 
non-thermal (magnetospheric) and thermal (polar-cap) 
components in the spectrum, potentially distinguishable 
from each other by phase-resolved spectroscopy of the 
source. 

The goal of our deeper XMM-Newton observation was 
to collect more photons for studying the pulsar spectrum. 
In addition, the high time resolution of the EPIC-pn de- 
tector could detect X-ray pulsations and potentially sep- 
arate its non-thermal and thermal components through 
phase-resolved spectroscopy. The observation and data 
analysis are described in Section 2, the relation between 
the detected X-ray and radio pulse shapes are described 
in Section 3, followed by a discussion of the results in 
Section 4. 

2. OBSERVATION AND DATA ANALYSIS 
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Fig. 1. — X-ray images (0.2-2.5 keV) of the field in the direction of 
the PSR J0108-1431. The pn image (~ 4' x2'), and the MOS1 and 
MOS2 images (each ~ 2.6' X 2.6') are shown on the top, bottom- 
left and bottom-right, respectively. North is up, East is to the left. 
The source extraction regions (radii of 12" for pn, 13" for MOS) 
are marked with blue circles, the background regions are marked 
with red circles. 
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Fig. 2. — Background light curve of the EPIC-pn detector in the 
energy range 10 to 12keV, binned to 100 s. The x-axis denotes 
the time counted from the obervation start, the y axis denotes the 
count rate in counts s —1 . The red solid line shows the considered 
exposure time range for the chosen GTIs for our timing analysis, 
requiring a count rate < 4 counts s . Similarly, the cyan line 
shows the chosen GTIs for our spectral analysis, requiring a count 
rate < 5 counts s —1 . Vertical lines indicate corresponding regions 
of flagged exposure times. 

Pulsar J0108-1431 was observed on 2011 June 15 
(MJD 55727) with XMM-Newton (obsid 0670750101) us- 
ing the European Photon Imaging Camera (EPIC) in 
Full-Frame mode with the Thin filter. The EPIC-pn 
(Striider et al. 2001) and EPIC-MOS1 and MOS2 cam- 
eras (Turner et al. 2001) observed the source for 126.7 
ks and 127.8 ks respectively. The EPIC data processing 
was done with the XMM-Newton Science Analysis Sys- 
tem (SAS) 11.0 2 , applying standard tasks. Our observa- 
tion was contaminated by soft proton flaring events in the 
background. We considered different selections of Good 
Time Intervals (GTIs) for optimal spectral and timing 
analysis, which are shown in Figure 2 and described in 
more detail below. 

We detected the X-ray counterpart of PSR J0108-1431 
at a position a = 01 h 08 m 08?25, 5 = -14°31'53'.'4 with 
EPIC-pn. This position is separated by 2'.'9 from 
the earlier (MJD 54136) Chandra detection at a = 
01 h 08 m 08^354, 5 = -14°31'50"38 (Pavlov et al. 2009) 
after accounting for the proper motion by Deller et al. 
(2009). The X-ray position differs by 3'.'0 from the ex- 
pected radio pulsar position using the proper motion and 
VLBI position: a = 01 h 08 m 08!347, 8 = -14°31'50'/187 
(MJD 54100) listed by Deller et al. (2009). The absolute 
astrometric 2a error of XMM-Newton is 4" 3 . Thus, the 
EPIC-pn X-ray position is consistent with the Chandra 
and VLBI radio position. 

2.1. Spectral Analysis 

For our spectral analysis of the X-ray counterpart of 
PSR J0108— 1431, we aimed for the highest possible 
signal-to-noise ratio, S/N = (N s - aN b )/(N s + a 2 N b )^, 
where N s and N b are the total numbers of counts 
extracted from the source and background regions of 

2 http://xmm.esac.esa.int/sas 

3 XMM calibration document: 

http: / /xmm2.esac.esa.int/docs/documents/C AL-TN-0018.pdf 
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Fig. 3. — Count rate spectra of the binned counts (minimal S/N> 
5) in the source region (black), background counts (blue circles), 
and background subtracted source counts (red) for EPIC-pn. 

areas A s and A b , respectively, and a = A s /A b . Given the 
strong background flaring in our observation (Figure 2), 
we tested different GTI filters derived from 100 s binned 
light-curves of events with energies above 10 and 12 
keV for pn and MOS detectors, respectively. For each 
GTI-filtcrcd event file, we applied the eregionanalyse 
task to optimize the aperture size for high S/N spectral 
extraction. Standard pattern filtering (< 4 for pn and 
< 12 for MOS) was enforced. We achieved the highest 
S/N for a GTI count rate cut-off of 5.0 counts s" 1 and 
0.5 counts s^ 1 for the pn and MOS light-curves, and for 
source aperture radii of 12" and 13" for pn and MOS, 
respectively. The net GTIs for pn, MOS1 and MOS2 
are 104.8 ks, 112.5 ks, and 114.6 ks, respectively. 

We used larger background regions for better statis- 
tics, with the radii of 42", 75" and 90" for pn, MOS1 
and MOS2, respectively. The extraction parameters are 
presented in Table 1. The source and background re- 
gions in the three processed EPIC images are shown in 
Figure 1. The background signal is comparable to the 
source signal at 1.2 to 2keV and then exceeds the source 
signal for energies > 3keV (Figure 3). Considering the 
substantial count rate errors for the high energies, we 
restrict our spectral fitting to events in the 0.2-2.5 keV 
range in order to better constrain the spectral fit param- 
eters. 

Considering all three EPIC cameras together, we 
obtained around 990 total X-ray counts from the source 
apertures, or around 705 source counts (see Table 1). 
Redistribution matrices and effective area files were 
generated using the usual SAS tasks rmf gen and arf gen. 
For the spectral fit, the SAS task specgroup was used to 
group the source counts of each spectra with a minimal 
S/N> 5 per bin. 

Using XSPEC 12. 6. 4 , we tested different spectral 
models (PL - powerlaw , BB - bbodyrad , BB+PL, 
BB+BB) applying % 2 statistic. For the photoelectric 
absorption in the interstellar medium (ISM), we used 
tbabs with the solar abundance table from Anders & 

4 http:/ /heasarc. gsfc.nasa.gov/docs/xanadu/xspcc 
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TABLE 1 

Spectral extraction parameters for pn, MOS 1/2 events. 



Extraction parameters 


pn 


MOS 1 


MOS 2 


Aperture radius 


12" 


13" 


13" 


Energy range (keV) 


0.2 - 2.5 


0.3 - 2.5 


0.3 - 2.5 


Total aperture counts 


682 


155 


152 


Source / total counts (%) 


69.9 


74.0 


74.4 


BG-cor. count rate (ks — 1 ) 


4.6 ±0.3 


1.0 ±0.1 


1.0 ±0.1 



Grevesse (1989) and the photoelectric cross-section table 
from Balucinska-Church & McCammon (1992) together 
with a new He cross-section based on Yan et al. (1998). 
We performed simultaneous fitting of the pn and MOS 
spectral data. 

The values of the fits with separate normalizations 
for each instrument agreed within the 90% confidence 
levels with the fit using the same normalization for all 
instruments. For simplicity, we give only the latter in 
the following . The spectral fit parameters were allowed 
to vary freely for the single component models, while 
for the combination of the PL with the BB there were 
not enough counts for constraining all the parameters 
sufficiently. Therefore we had to freeze a parameter (see 
below). 



An absorbed 

p _ o 1+0.5 

N H = 



(5.51*1) x 10 



PL 

and 

20, 



model with photon index 
Hydrogen column density 
cm~ 2 provides a good fit of 
the data with xl = 1-1; see Figure 4 for the X-ray spec- 
tral fit and Figure 5 for confidence contours for the three 
model parameters. Table 2 lists our best spectral fitting 
parameters. The confidence levels for each parameter 
were obtained using XSPEC's error and steppar 
commands. The non-thermal luminosity is estimated to 



hp r PL 

De ^0.2-2.5 keV 



lira ^0.2-2.5 kcV 



x 10 29 erg 



s . The luminosity uncertainties include both model 
and distance uncertainties. 

In contrast to the PL fit, a pure thermal model does 
not describe the data well. The best fit of the absorbed 
BB model has large residuals (xt = 2.3) that rule out 
the possibility of the emission being entirely BB-like. In 
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Fig. 4. — Absorbed PL fit and its residuals in units of sigmas for 
the pn (black), MOS1 (red) and MOS2 (green) data. 




Fig. 5.— 68%, 90% & 99% confidence contours in the T - 
Nh plane (top) and in the T — PL normalization plane. Nh 
is units of 10 21 cm - 2 , the PL normalization J\f— a is in units of 
10 — 6 photons cm -2 s — 1 keV - 1 measured at IkeV. Plotted are also 
the contours of constant unabsorbed flux, F uuahB in units of 10 — 14 
erg cm" 2 s _1 for the 0.2 - 2.5 keV band (bottom). 



TABLE 2 

Best fit parametric values with 90% confidence limits. 



Parameters 



Best fit values 



Absorbed PL fit 



N H (10 20 cm" 2 ) 

r 

PL norm (10 — 6 photons cm -2 s — 1 keV - 1 at IkeV) 
X 2 /d.o.f. 

F^ s „ (10~ 15 ergon" 2 s" 1 ) 



0.2-2.5 koV 
punabs 
^0.2-2.5 keV 



(10 



-14 —2 

erg cm i 



- 1 ) 



r + l.b 

'-2.2 
I +0.5 
-0.2 
ii+0.4 
-0.4 

1.1/27 
9 4+0-9 

M '^-1.0 



5.5" 1 
3.lj 
2.8"* 



2.3 



+1.2 

-0.8 



Absorbed PL+BB fit 



N H (10 



20 



PL norm (10~ 6 photons cm -2 s _1 keV" 1 at IkeV) 
kT (keV) 
/?bb a (m) 
Xl/d.o.f. 

F* b o s „ (10- 15 ergcm- 2 s- 1 ) 



— 

-2.30 
2.0 Fixed 
7+0-4 
-0.5 
,+0.03 
-0.01 
3 + 16 
-9 



2.3: 



1.7Z 
0.11 H 

43"* 



1 0.2-2.5 keV 
punabs (1(1 
- r 0.2-2.5kcV l ±u 



- 1 ) 



1.1/26 

9 7+ 10 
M -'-1.0 

1.3 " 



+0.5 

0.3 



a The radius of the blackbody emission was obtained from the 
normalization using a distance of d = 210 pc (Verbiest et al. 2012). 
Its error is the propagated normalization error only. If one also 
considers the distance error the correponding values are i?t,b = 



43: 



principle, one could try to fit the spectrum with neutron 
star atmosphere (NSA) models (Zavlin et al. 1996; 
Pavlov et al. 1995). However, the available atmosphere 
models in XSPEC were calculated cither for B = 
or strong magnetic fields (B = 10 12 G and higher). 
For the magnetic field strength of PSR J0108-1431, 
B = 2.5 x 10 11 G, the (redshifted) cyclotron energy is 
~ 2 — 3keV. This is too close to the observed photon 
energy range and has a strong effect on the atmophere 
spectrum, which is not considered in the weak or strong 
magnetic field models in XSPEC. Therefore, the NSA 
models are not applicable for PSR J0108-1431, and we 
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Fig. 6. — Top: Absorbed PL+BB fit and the residuals in units 
of sigma deviations, pn (black), MOS1 (red) and MOS2 (green). 
Bottom: Underlying individual PL and BB spectral model compo- 
nents. 

have to stick to the simplistic BB models. 

We also checked a two-component BB+BB model 
with photoelectric absorption. Such a model could 
describe the thermal emission from a nonuniformly 
heated neutron star surface. The best fit for this model 
also has residuals too large to be acceptable (xt = 2.0). 

A combination of non-thermal magnetospheric emis- 
sion and thermal emission from polar caps can be 
modeled by a PL+BB model. Because of the noise in 
the data, we had to freeze one fitting parameter. We 
froze the photon index to T = 2, motivated by the 
typical photon indices of younger pulsars, r ~ 1 — 2 
(Li et al. 2008). The fit is acceptable (xl = 1.1), its 
BB temperature is kT = 0.11^°;^ keV. We use the 
bbodyrad model in XSPEC, whose normalization factor 

gives the projected area Aj_ = 5700^2400 ^210 m2 an( l 
the corresponding radius of an emitting equivalent 
sphere, i?bb = (A^/n) 1 / 2 = 43 ig 6 d2io m, where <i 2 io is 
the distance divided by 210 pc (the errors include only 
the propagated normalization error), see Table 2 for 
additional fit results. Since i?NS 3> i?bt>, the emission 
likely comes from part of the surface like hot polar 
caps. Figure 6 shows the spectral fit and its residuals, 
while the 68%, 90%, and 99% confidence contours in 
the kT-Nn and kT-A± planes are plotted in Figure 7. 
As mentioned above, there are no available atmosphere 
models for the magnetic field strength of PSRJ0108- 
1431. Therefore, we do no fit a two-component PL+NSA 




0.12 
kT (keV) 

Fig. 7.— 68%, 90% & 99% confidence contours in the kT - N H 
plane (top) and kT - Aj_ plane (bottom) for the absorbed PL+BB 
model. Nh is in units of 10 21 cm -2 . Lines of constant bolometric 
luminosities in units of 10 28 ergs -1 are overplotted in the kT - Aj_ 
plane. 

model to the data. 
The luminosity for both components together is esti- 



mated to be L^ L 5kcV 



punabs 

— ^0.2-2.5 koV 



And 2 = 6.8±|;f 



10 28 erg s 1 . The luminosity of the non-thermal compo- 



nent is Lj^_ 2 . 5kcV = 3.7±|;? x 10 28 erg s -1 



The derived 

values of the temperature and area correspond to the 
bolometric thermal luminosity of an equivalent sphere 
£bb = 4A ±a T 4 = 3.7±|;f x 10 28 erg s _1 . The lumi- 
nosity errors include the propagated distance error and 
the corresponding propagated fit errors in flux, temper- 
ature and emission area. Lines of constant bolometric 
luminosities are overplotted in the kT - A± plane (Fig- 
ure 7, bottom). The above listed errors of Lboi encom- 
pass nearly the whole 3a contours in the kT - A± plane 
because we additionally considered the distance error in 
our conservative error estimate. 

2.2. Timing Analysis 

Hobbs et al. (2004) reported the radio ephermerides 
of PSR J0108-1431 to be v = 1.23829100810(3) Hz and 
v = 0.11813(18) x 10- 15 Hzs _1 at MJD 50889.0. At the 
start time of our X-ray observations, MJD 55727.2, the 
expected frequency change, — 0.049 //Hz, is negligible for 
X-ray timing. The EPIC-pn detector with a nominal 
frame time of 73.4 ms in the full frame mode is well 
suited for the timing analysis of PSR J0108— 1431, while 
the time resolution of the MOS detectors in full frame 
mode, 2.6 s, does not allow such an analysis. 

Following up on SAS warnings during the data process- 
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ing, we checked the EPIC-pn CCD with the target on it 
for time jumps in the data. We set the SAS environment 
parameter SAS.JUMP .TOLERANCE 5 to a value of 33 
to account for the known temperature effect on the ac- 
tual frame time of the EPIC-pn detector (Freyberg et al. 
2005, Freyberg et al. 2012 6 ). In addition, we excluded 
the end of the observation (> 119.28 ks after start). At 
that time, a bright X-ray background flare caused the in- 
strument to switch to the counting mode 7 which resulted 
in finding of false time jumps by the SAS. Our final data 
for the timing analysis were free of apparent time jumps. 
All times-of-arrival of the X-ray photons were corrected 
to the solar barycentric system using the standard task 
barycen. 

We used the Z 2 test, the sum of powers of first n 
harmonics, (e.g., Buccheri et al. 1983) to search for 
pulsations in these data. For our timing analysis, we 
checked 7 different GTI screenings, 11 different energy 
regions, and 5 different extraction regions in order to 
maximize the Z\. The GTI screening with pn count rate 
< 4 counts s -1 for a 100 s binned background light curve 
(10-12 keV), the extraction radius of 8", and the energy 
range of 0.15-2 keV were found to be the optimal choices 
among the tested variants. At the chosen GTI filter, 
there is only one small (200 s) time gap in the last fifth 
of the exposure; see Figure 2. This filtering provided 
507 events during a time span of Tspan = 119.1 ks. 

PSR J0108— 1431 has been monitored for 16 years, and 
no glitches have been detected (Espinoza et al. 2011). In- 
vestigating possible timing irregularities for 366 pulsars, 
Hobbs et al. (2010) found PSR J0108— 1431 to be a very 
stable pulsar with low timing noise like other pulsars with 
similarly low v. Therefore, we could search for X-ray 
pulsations at the radio pulsation frequency. However, 
to check whether the above-mentioned time jump cor- 
rection worked properly, we searched a frequency range 
of 1.2382-1. 2384Hz with a sampling of 0.1 //Hz, thus 
oversampling the expected Z 2 peak width, ~ Tj,"^ = 
8.4 fiHz, by a factor of 80. A peak Zf max = 48.0 is 
found at vq — 1.2382908 Hz ±0.7/iHz, corresponding to 
a period of 0.8075647s ±0.5/j,s. The frequency uncer- 
tainty was derived similarly to Chang et al. (2012) as 
Sv = 0.55T sp i n (Z 1 2 max )- 1 / 2 . Within errors the X-ray 
pulse frequency agrees with the radio pulse frequency 
very well. The probability to find Zf max = 48.0 by 

chance is p = exp(— Z\ max /2) — 4 x 10 , which corre- 
sponds to the 6.6tr confidence level. 

To explore the harmonic content of the X-ray pulsa- 
tions and refine the significance estimate, we performed 
the Z 2 test accounting for up to n = 7 harmonics. The 
peak Z 2 values (and the respective significances) are 

61.5 (6.7cr), Z\ = 66.4 (6.7a), 



zl 



57.9 (6.8a), Z\ 



5 SAS-JUMPJTOLERANCE is given in units of 20.48 ms. Only 
deviations of the measured actual frame time from the nominal 
frame time larger than the SAS.JUMP_TOLERANCE are identi- 
fied as time jumps by the SAS. 

6 See also: http://www2.le.ac.uk/departments/physics/research/ 
src/Missions / xmm-newton/technical /leicester-2012-03/freybcrg- 
cal-2012.pdf 

7 XMM-Newton Users Handbook, section 3.3.2 - Science modes 
of the EPIC cameras 




1.0 

pulse phase 

Fig. 8. — Folded light curves as histograms and reference phase 
averaged (RPA) pulse profiles (see text). For a better overview, 
two full phase cycles are shown. 

Upper panel: The upper If -bin histogram represents the full 
considered energy range, 0.15-2 keV, the corresponding RPA pulse 
profile is shown in blue, its average background level is indicated 
with a blue dashed line. Similarly, for the energy range 0.15—1 keV, 
the red RPA pulse profile corresponds to the lower histogram, and 
the average background level for this energy range is indicated with 
a red dashed line. Lower panel: The 6-bin histogram and the 
RPA pulse profile represent the folded light curve in the energy 
range 1-2 keV. 



Zl = 75.3 (6.9ct), Zl = 75.5 (6.6a), Z 2 = 80.0 (6.4cr) at 
frequencies consistent within errors with the one found 
from the Z\ test. For the application of the H-tcst 
(de Jager et al. 1989), H = ma.x{H m }, we found the 
H m = Z? n - 4m + 4 to be 48.0, 53.9, 53.5, 54.4, 59.3, 
55.5, 56.0 for the first to seventh harmonic, respec- 
tively, i.e. the pulsation have statistically significant 
contributions from 5 leading harmonics of the principal 
frequency. 8 Thus, the pulsations of PSR J0108-1431 are 
unambiguously detected in X-rays with a siginificance 
of6.9cr. 

In order to visualize the pulse shape, we obtained a 
folded light curve in form of a histogram for the energy 
range 0.15-2 keV (507 counts), as well as for the splittcd 
energy ranges 0.15-1 keV (416 counts) and 1-2 keV (92 
counts). The frequency Vq = 1. 2382908Hz, was used for 
the folding. We applied the so-called Scott's rule for 
setting the upper bound for the histogram bin width by 
Terrell & Scott (1985). These authors concluded that 
their rule to avoid oversmoothing gives nearly optimal 
results for a variety of smooth probability densities. 
However, choosing the optimal histogram bin size 
remains a matter of debate in statistical data analysis 
(for a review, see e.g., Scott 1992). According to Scott's 
rule the number of bins must be > (2./V) 1 / 3 , where N 
is the number of events. We selected JVhist = H bins 

8 We checked H m up to 20 harmonics as recommended by de 
Jager et al. (1989) with the same result. 
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for the energy ranges 0.15-2 keV and 0.15-1 keV, and 
•Whist = 6 for the energy range 1-2 keV. 

The number of counts in the histogram bins depends 
on the reference phase, which was chosen arbitrarily. 
The histogram can look quite different for another 
choice of reference phase, especially if there are few 
bins. To obtain the folded light curve independent of 
the reference phase choice, we averaged the histogram 
over the reference phase 9 , varying the latter within one 
histogram bin; see the Appendix A for a short descrip- 
tion of the used algorithm. In the following, we call the 
obtained new histogram the reference phase averaged 
(RPA) pulse profile. The histograms of the folded light 
curve together with the respective RPA pulse profiles 
are shown in Figure 8. The pulse shape for the full 
energy band, 0.15-2 keV, and the soft energy band 
appear to be asymmetric, with a slower rise and steeper 
decay. The two RPA pulse profiles have their maxima at 
similar phases. The 1-2 keV RPA pulse profile appears 
more sinusodial and slightly shifted to smaller phases 
compared to the low-energy and broad bands. However, 
the pulse shape in the energy range 1-2 keV has poorer 
statistics, and direct comparisons must be regarded with 
caution. Using the RPA pulse profiles as probability 
distribution templates and the respective number of 
events as sample sizes, we carried out Monte Carlo 
simulations of the profiles to estimate the uncertainties 
of the phase positions of the respective RPA pulse 
profile maxima. The maxima of the RPA pulse profiles 
in Figure 8 are at 4> = 0.96 ± 0.03, 0.97 ± 0.03, and 
0.86 ±0.06 for the energy ranges 0.15-2 keV, 0.15-1 keV, 
and 1-2 keV, respectively (see Figure 8). 

For the full energy band, 0.15-2 keV, for which we 
have the best statistics, we also use another visualiza- 
tion method. We apply Fourier decomposition F m (<j>) up 
to an harmonic m to describe the pulse shape. 



F m {<t>) 



where 



N 



hist 



1 + 2 a,k cos(2irk(f)) + bk sm(2%k(j)) 



fe=i 



(1) 



1 N 1 N 

afe = — ^cos(27rfc^), b k = — 2jsm(27rfci>ij) (2) 

»=i i=i 

are one half of the empirical Fourier coefficients, <j> is 
the phase, N is the total number of events in the 
folded lightcurve, U are the event times, and v = vq = 
1.2382908 Hz is the pulsation frequency. The scaling fac- 
tor iV/iVhist allows one to compare the Fourier series plot 
with the iVhist-bin histogram. Here, we apply a scaling 
factor N/N hist = 507/11 « 46 to plot the Fourier se- 
ries for m =3, 4 and 5 in Figure 9. While including five 
harmonics introduces 'wiggles' of the same order as the 
errors in a light curve bin, Fourier series with coefficients 
including up to three or four harmonics appear to model 
well the folded light curve, in particular its asymmetric 
pulse shape. 




pulse phose 

Fig. 9. — Folded light curve histogram for the full considered 
energy range 0.15— 2 keV together with the Fourier series including 
up to three harmonics (dark blue), four harmonics (cyan) and five 
harmonics (yellow). For comparison, the RPA pulse profile from 
Figure 8 is also plotted as red dashed curve. 



From the histograms, we calculate the empirical pulsed 
fractions as the ratio of the number of counts above the 
minimum to the total number of counts. For the energy 
ranges 0.15-2 keV, 0.15-1 keV, and 1-2 keV the pulsed 
fractions are f p = 39 %±6 %, 33 %±7 %, and 80 %±15 %, 
and the background-corrected, instrinsic pulsed fractions 
are f p nt = 47%±7%, 38%±8%, and 100 ±° 1S % respec- 
tively. We calculated the errors of the pulsed fraction as 
Sf p = (2/jV)5 j following the estimate for broad pulses by 
de Jager (1994) 10 . 

3. COMPLEMENTARY RADIO DATA 

It is interesting to compare the radio pulse with the 
X-ray pulse, in particular the phase difference between 
them. We use the monitoring radio data described by 
Weltevrede et al. (2010) for our comparison of the ra- 
dio and the X-ray profiles. Observations at 1.4 GHz were 
obtained from December 2009 to July 2012 using the 
64-m radio telescope at Parkes, NSW, Australia. The 
data were calibrated and times-of-arrival (TOAs) pro- 
duced using psrehive routines (Hotan et al. 2004) and 
the tempo2 package (Hobbs et al. 2006). An accurate 
value of the radio dispersion measure (DM) is neces- 
sary to correct for the dispersion delay when comparing 
the relative phases of the 1.4 GHz radio and high-energy 
pulse profiles. The measured DM for PSR J0108-1431 is 
2.38 ± 0.19 cixr 3 pc (Hobbs et al. 2004). This translates 
into an DM-induced uncertainty of ~ 402 fis or 5 x 10~ 4 
in the phase of the radio pulse. 

For the comparison with the radio pulse shape, we con- 
verted the TOAs of the X-ray photons to the solar sys- 
tem barycenter using the standard SAS task barycen 
with the same coordinates as in the radio tempo2 anal- 
ysis and the same DE405 solar-system ephemeris (Stan- 
dish 1998). We used the barycentric TOAs of five X -ray 
events in tempo2 to derive their phases relative to the 
radio profile. Knowing the X-ray phases of these events 
as well, we determined the offset between the radio and 
X-ray phase reference points to be 0.898 ± 0.003. Cor- 
recting for this shift between the reference systems, Fig- 
ure 10 shows the X-ray and radio profiles together in the 
same phase range. We confirm the unusual high linear 



9 A similar phase averaging was applied by Zavlin et al. (2002) 
(their Figure 5). 



10 As the power in the fundamental harmonic strongly exceeds 
the powers in higher harmonics for PSR J0108— 1431, this uncer- 
tainty estimate is applicable, albeit approximate. 
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Fig. 10. — The folded X-ray light curve histogram (energy range 
0.15-2 keV) is plotted together with the corresponding RPA pulse 
profile (both in black). The red curve shows the total intensity of 
the 1.4 GHz radio pulse in mjy; the blue curve shows the linearly 
polarized part of the radio pulse. The reference phases of the X-ray 
and radio profiles were aligned as described in the text. 



polarization of the radio pulse reported by Weltevrede 
et al. (2010), obtaining a value of 77% for the fractional 
linear polarization at the pulse peak. Using the peak 
positions of the RPA pulse profiles from Section 2.2, the 
radio pulse leads the X-ray pulse peak by 0.06 ±0.03 and 
0.08 ± 0.03 in the energy ranges 0.15-2 keV and 0.15- 
1 keV, respectively. In contrast, the X-ray pulse at 1- 
2keV leads the radio pulse by 0.03 ± 0.06. The quoted 
errors take into account the uncertainties in the reference 
phase shift, the DM-induced phase uncertainty, the ab- 
solute XMM-Newton EPIC-pn timing accuracy of 48 fj,s 
(Martin-Carrillo et al. 2012), and the dominating error 
of the maximum position in the RPA pulse profile (Sec- 
tion 2.2). 

4. DISCUSSION 

4.1. The X-ray spectrum and luminosity 

PSR J0108— 1431 is the oldest among the non-recycled 
ordinary pulsars detected in X-rays, with the spin-down 
age a factor of 4 larger than that of PSRB1451— 68, the 
previous record holder. It also has the lowest E among 
the same sample of X-ray detected non-recycled pulsars. 
Our new XMM-Newton observation provided a factor of 
20 more source counts enabling better characterization 
of the pulsar's spectrum compared to the previous 
Chandra observation. However, the strong flaring 
and high background at energies above ~ 2 keV have 
undermined energy-resolved timing and phase-resolved 
spectral analysis of the pulsar emission. 

Formally, a simple absorbed PL model and an ab- 
sorbed PL+BB model describe the data equally well. 
The photon index V — 3.1 to' 2 ( see a ls° Figure 5) of 
the PL fit is larger (i.e., the spectrum is steeper) than 
T ~ 1-2, typical for young pulsars (Li et al. 2008). As 
mentioned above, for the BB+PL fit we had to freeze 
one parameter because of the small number of counts 
and strong noise. We chose to fix the photon index 
at r = 2, assuming similar magnetospheric emission 
characteristics for young and old pulsars. 

The PL fit falls in line with the trend observed in 



other old pulsars, for which the PL fits suggest too 
large Nh values in conjunction with larger photon 
indices T. For the line of sight of PSR J0108-1431, 
the LAB Survey of Galactic neutral hydrogen reports 
Nhi = 2.1 x 10 20 cm~ 2 in this direction (Kalberla 
et al. 2005), the Dickey & Lockman (1990) neutral 
hydrogen survey reports Nhi = 1-8 x 10 20 cm~ 2 . The 
Nh = 5. 5^2 2 x 10 20 cm~ 2 from our PL fit is above 
these total Galactic values. We use the Lutz-Kelker 
corrected, parallactic distance d = 210^5° P c by Ver- 
biest et al. (2012) to estimate the expected Nh value 
applying the 'analytical' 3D extinction model described 
by Posselt et al. (2007). For the close (< 250 pc) solar 
neighbourhood, this model is based on the 3D NaD 
absorption line mapping by Lallement et al. (2003) 
and has a resolution of ~ 25 pc; at larger distances an 
analytical model is used for the extinction (see Posselt 
et al. 2007, 2008 for more details). Considering the 
errors of the distance, we derived an expected Nh is 
in the range (0.3-0.8) x 10 20 cm~ 2 . Assuming 10 H 
atoms per electron, we derive a similarly low expected 
Nfl M — 0.7 x 10 20 cm~ 2 from the pulsar dispersion 
measure, DM = 2.38 pccm~ 3 (Hobbs et al. 2004). Thus, 
the absorption is overestimated in the simple PL fit (see 
Figure 5) while the 90% confidence range of the Nh 
parameter in the BB+PL fit includes the expected Nh 
value. 

The temperature obtained in the BB+PL fit, 
T = 1.3 MK, is a reasonable estimate for the expected 
heated polar cap region. The effective projected 
emitting area is A± = 5700 m 2 . This is a factor of 
~ 14 smaller than the conventional polar cap area 
A pc = 2ir 2 R^ s /{cP) « 82000m 2 , using R NS = 10km. 
PSR J0108— 1431 is very similar to other old pulsars 
in this respect (Posselt et al. 2012; Pavlov et al. 2009; 
Misanovic et al. 2008; Kargaltsev et al. 2006; Zhang 
et al. 2005). We should note, however, that the values 
of the temperature and, particularly, the projected area 
are rather uncertain because of their strong correlation 
(Figure 7). 

As mentioned in the Introduction, one has to take 
into account that thermal emission from the surface 
layers of a neutron star can differ substantially from the 
BB emission (e.g., Pavlov et al. 1995). In particular, 
if the emission emerges from a light-element (H or 
He) atmosphere, fitting with neutron star atmosphere 
models may yield the effective temperature a factor of 
2 lower than Tm, and the projected emitting area a 
factor of 10-100 larger than A±, while the bolometric 
luminosity does not change substantially. The magnetic 
field can have a strong impact on the observed emission 
of a neutron star atmosphere. The current neutron star 
atmosphere models available in XSPEC are not applica- 
ble to PSR J0108-1431 because the electron cyclotron 
line, caused by its magnetic field, B = 2.5 x 10 11 G, is 
in the observable energy range, but it is not included 
in the XSPEC models. Therefore, we do not discuss 
atmosphere models for PSR J0108-1431. 

Overall, there are now several old pulsars whose 
spectra could be described by the PL model, but the 
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Fig. 11. — X-ray luminosities and upper limits of eleven old pul- 
sars versus their spin-down energies E. In cases of sufficient num- 
bers of counts, the diamonds and asterisks show the 1-10 keV 
non-thermal luminosities and bolometric thermal luminosities, re- 
spectively; otherwise the luminosities were obtained from PL fits. 
The arrows mark upper limits derived from X-ray non-detections. 
This figure is an update of the Figure 4 presented by Posselt et al. 
(2012) using distances corrected for the Lutz-Kelker bias by Ver- 
biest et al. (2012). In case of PSRB1451-68, error propagation 
leads to a lower error of the bolometric thermal luminosity as large 
as the best-fit value. Therefore, we plot only an upper limit for the 
bolometric thermal luminosity of PSRB1451-68. 

properties of these fits - in particular, the large T 
and the much higher than expected Nh - suggest a 
more complicated model. Adding a thermal (e.g., BB) 
component to the spectral model allows lower N H values 
and, correspondingly, smaller T, similar to those of 
younger pulsars. 

For comparison with other pulsars, we calculated the 
non-thermal (PL) luminosity in the 1-10 keV energy 
= 2.0 



range: I^-iokeV 



x 



10 28 ergs" 



PL 



-lOkoV 



= 3.3 



+i. 

l ° 9 x lO^ergs- 1 



1 for the PL 
for the PL 



fit, and L\ 

component in the BB+PL fit 11 . The higher photon 
index of the PL fit, obtained from the fitting at a lower 
energy range, is the reason why the L PL is smaller in the 
case of the PL fit than in the case of the PL+BB fit. The 
errors were calculated from the 90% confidence errors of 
the unabsorbed flux and the distance uncertainty of the 
Lutz-Kelker corrected distance as listed by Verbiest et al. 
(2012). These luminosities translate into non-thermal X- 



ray efficiencies of ryf 



PL 



-lOkcV 



-^l-lukcy/^ 



0.003 d 2 . 



and - 0.006 d 2 2W 
fit, respectively. 



rbb 



— q 7+5.8 

— d-l_ 2 .7 



for the PL fit and the PL+BB 
Using the bolometric luminosity 
10 28 erg s _1 (Section 2.1), one can 
also estimate the thermal polar cap heating efficiency 
from the PL+BB fit, rfe h c = L^JE ~ 0.006 d| 10 , 
which gives the fraction of spindown power heating the 
polar caps. Harding & Muslimov (2001) predicted that 
the expected polar cap heating efficiency for ordinary 
pulsars for ages r < 10 7 years grows with age and period 
(e.g., rjp C ~ 0.01 for r = 1.5 x 10 7 years, P = 0.2 s; 
see their Figure 7). For pulsars with r > 10 7 years 

11 The 1-10 keV PL component luminosity and its errors were 
estimated for the fixed T = 2. 



they cautioned that the pulsar cannot produce enough 
electron-positron pairs to fully screen the electric field. 
Thus, the fraction of the returning positrons that heat 
the polar cap decreases, and the heating efficiency 
will drop. Comparing the estimated polar cap heating 
efficiency with the predictions by Harding & Muslimov 
(2001), PSRJ0108-1431 has a lower efficiency than 
one would expect for a pulsar having the same period 
(0.8 s) at an age of r = 10 7 years. This could indicate 
that the returning positron fraction is indeed smaller for 
PSR J0108-1431 than for younger pulsars. 

In Figure 11, we updated the X-ray luminosities ver- 
sus spin-down power plot by Posselt et al. (2012) (their 
Figure 4 12 ). The pulsar again shows properties compara- 
ble to those of other old pulsars. In particular, its non- 
thermal X-ray efficiency is higher than those of young 
pulsars, most of which have rjfz w]ic y < 10~ 3 . Thus, the 
observation that older pulsars seem to radiate in X-rays 
more efficiently than younger ones (Zharikov et al. 2004; 
Kargaltsev et al. 2006) is reinforced. 

4.2. The X-ray pulsations of PSR J0108-1431 

The timing analysis of PSR J0108-1431 establishes, 
for the first time, unambiguous X-ray pulsations from 
this pulsar. More than 2 harmonics are required to 
explain the asymmetric pulse profile. Such asymmetric 
pulse profiles can have several explanations, different for 
nonthermal and thermal emission. 

Asymmetric pulse profiles can be easily produced by 
nonthermal emission in the outer gap model (Romani 
& Yadigaroglu 1995), especially if the special relativity 
effects, such as aberration and retardation, are taken 
into account (e.g., Dyks & Rudak 2003 and references 
therein). A comparison of the detected pulse shapes 
with those predicted by the models of magnetospheric 
high-energy emission could be useful in establishing the 
geometry of the emitting region, especially if additional 
information about directions of the pulsar's magnetic and 
spin axes is available from, e.g., radio-polarimetry. For 
PSR J0108-1431, we do not know if the bulk of the ob- 
served X-ray emission is produced in the pulsar's mag- 
netosphere (as it is implied by the PL spectral model), 
or whether only a high-energy tail originates from the 
magnetosphere (as we assumed in the PL+BB model). 
A narrow, sharp pulse would unambiguously point to 
strongly beamed magnetospheric emission; however, the 
pulse profile appears to be rather broad (Figure 8) even 
at higher energies (1-2 keV). Thus, we can neither rule 
out nor confirm purely non-thermal emission as the cause 
for the asymmetric pulse profile. 

The asymmetric pulse profile can also be produced 
by thermal emission if the distributions of temperature 
and/or magnetic field over the neutron star surface are 
not axially symmetric. For instance, an asymmetric 
pulse with a shape similar to the observed one could 
be produced by a polar cap whose trailing edge is 
hotter (hence brighter) than the leading edge, as it 

12 Note that Posselt et al. (2012) used distances derived from 
the Lutz-Kelker corrected parallaxes (Verbiest et al. 2010), which 
differ from the Lutz-Kelker corrected distances. See Verbiest et al. 
(2012) for details about the differences. 
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was suggested by McGowan et al. (2003) to explain 
a similar pulse shape for the young (r = 30 kyr) 
PSR J0538+2817. Another possibility to explain the 
asymmetric pulse of thermal radiation is provided by the 
anisotropy in the local intensity of radiation emerging 
from a neutron star atmosphere. In contrast to the BB 
radiation, atmospheric radiation in a strong magnetic 
field is highly anisotropic, with a 'pencil component' 
along the magnetic field, and 'fan component' transverse 
to the magnetic field (Pavlov et al. 1994). If the mag- 
netic field is not perpendicular to the polar cap surface 
(which can occur if the field is essentially nondipolar), 
then a great variety of pulse shapes can be produced, 
depending on the magnetic field inclination and the 
angle between the spin axis and line of sight. The 
strong anisotropy of the atmospheric emission could also 
explain the relatively high pulsed fraction of the thermal 
emission from PSR J0108- 1431. The pulsed fraction 
fp nt — 38 %±8 % in the energy range 0.15-1 keV appears 
to be too high for locally isotropic BB emission, whose 
pulsations would be strongly smeared by the bending 
of photon trajectories in the gravitational field of the 
neutron star (Zavlin et al. 1995). Thus, the properties 
of the pulsations detected in PSR J0108-1431 do not 
contradict to the presence of a thermal component in its 
emission. To firmly prove the presence of this component 
and infer the properties of the polar cap(s) from the 
comparison with the models, a refined energy-resolved 
timing analysis is needed, which would require "cleaner" 
data. 

We found a possible phase shift A<f) = 0.06 ± 0.03 be- 
tween the 1.4 GHz and 0.15-2 keV pulse peaks. In prin- 
ciple, it might be due to different emission heights, but 
the uncertainty of A</> is too large to warrant such an 
investigation. It is interesting to note, however, that the 
1-2 keV RPA pulse profile peak appears to be closer to 
the radio peak than the 0.15-1 keV RPA pulse profile 
peak. The 1-2 keV pulse seems to be slightly shifted 
by Acf> ~ 0.1 from the 0.15-1 keV pulse (considering the 
maxima of the RPA pulse profiles). In addition, as is 
apparent from Figure 8, the pulse profile at 1-2 keV ap- 
pears to be sinusodial while the soft energy pulse pro- 
file (0.15-1 keV) is asymmetric. Phase shifts and differ- 



ent pulse shapes would support the hypothesis of two 
components - an (anisotropic) thermal component and a 
non-thermal component. However, the high background 
and the associated large errors prohibit any firm conclu- 
sion. It remains to be confirmed, whether or not the 
two-component model is a viable interpretation, in par- 
ticular, if a separate non-thermal X-ray component is 
formed in close vicinity to the radio emission. 

5. SUMMARY 

We detected for the first time X-ray pulsations of 
PSR J0108— 1431 at the pulse frequency expected from 
radio pulse timing. The pulse shape is rather asymmet- 
ric, requiring up to 5 harmonics to describe it. The peaks 
of the radio and the X-ray pulses are close to each other 
in phase. The X-ray spectrum and the high non-thermal 
X-ray efficiency of PSR J0108— 1431 are comparable to 
other old pulsars. In particular, while the spectrum can 
formally be well described by a PL fit, the expected 
smaller photon index and lower Nh are better accounted 
for in a PL+BB model. We were not able to investigate 
phase-resolved spectra, because of the high background 
in the strongly flare-impaired observation. 
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APPENDIX 

OBTAINING THE REFERENCE PHASE AVERAGED PULSE PROFILE 

In the following, we briefly describe the sliding average bin algorithm used to average over reference phases within 
one histogram bin, producing the reference phase averaged (RPA) pulse profile. First, we create k = 0, 1, . . . , (A su b — 1) 
folded light curve histograms with Ahist bins by shifting the reference time in the time series by Atk = fcP/(Ahi s t x A su b), 
where P is the period of the pulsations, and N su b > 1 is an integer, which is equivalent to shifting the reference phase 
by A(j}}~ = Atk/P- As usual, the bin heights of the histograms are the number of photon counts in the respective bins. 
Second, we divide the whole phase interval into n = 1, 2, . . . A su b x Ahist sub-intervals with bin heights /i^, i.e., there 
are A su b sub-intervals in each original bin having all the same bin heights. Finally, we average the bin heights of the 
A su b histograms in each of the sub-intervals: 

1 iV sub -l 

k * FA = W- E k n ( A1 ) 
sub fc=0 

For our histograms and RPA pulse profiles in Figure 8, we used A su b = 100, and, as mentioned in Section 2.2, 
we used Ahist = 11 bins for the energy ranges 0.15 to 2 keV and 0.15-1 keV, and Ahi s t = 6 for the energy range 1-2 keV. 
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